home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / dht / dht.c next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  4.9 KB  |  226 lines
  1. /* dht/dht.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman
  21.  */
  22. #include <config.h>
  23. #include <stdlib.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_math.h>
  26. #include <gsl/gsl_sf_bessel.h>
  27. #include <gsl/gsl_dht.h>
  28.  
  29.  
  30. gsl_dht *
  31. gsl_dht_alloc (size_t size)
  32. {
  33.   gsl_dht * t;
  34.  
  35.   if(size == 0) {
  36.     GSL_ERROR_VAL("size == 0", GSL_EDOM, 0);
  37.   }
  38.  
  39.   t = (gsl_dht *)malloc(sizeof(gsl_dht));
  40.  
  41.   if(t == 0) {
  42.     GSL_ERROR_VAL("out of memory", GSL_ENOMEM, 0);
  43.   }
  44.  
  45.   t->size = size;
  46.  
  47.   t->xmax = -1.0; /* Make it clear that this needs to be calculated. */
  48.   t->nu   = -1.0; 
  49.  
  50.   t->j = (double *)malloc((size+2)*sizeof(double));
  51.  
  52.   if(t->j == 0) {
  53.     free(t);
  54.     GSL_ERROR_VAL("could not allocate memory for j", GSL_ENOMEM, 0);
  55.   }
  56.  
  57.   t->Jjj = (double *)malloc(size*(size+1)/2 * sizeof(double));
  58.  
  59.   if(t->Jjj == 0) {
  60.     free(t->j);
  61.     free(t);
  62.     GSL_ERROR_VAL("could not allocate memory for Jjj", GSL_ENOMEM, 0);
  63.   }
  64.  
  65.   t->J2 = (double *)malloc((size+1)*sizeof(double));
  66.  
  67.   if(t->J2 == 0) {
  68.     free(t->Jjj);
  69.     free(t->j);
  70.     free(t);
  71.     GSL_ERROR_VAL("could not allocate memory for J2", GSL_ENOMEM, 0);
  72.   }
  73.  
  74.   return t;
  75. }
  76.  
  77. /* Handle internal calculation of Bessel zeros. */
  78. static int
  79. dht_bessel_zeros(gsl_dht * t)
  80. {
  81.   unsigned int s;
  82.   gsl_sf_result z;
  83.   int stat_z = 0;
  84.   t->j[0] = 0.0;
  85.   for(s=1; s < t->size + 2; s++) {
  86.     stat_z += gsl_sf_bessel_zero_Jnu_e(t->nu, s, &z);
  87.     t->j[s] = z.val;
  88.   }
  89.   if(stat_z != 0) {
  90.     GSL_ERROR("could not compute bessel zeroes", GSL_EFAILED);
  91.   }
  92.   else {
  93.     return GSL_SUCCESS;
  94.   }
  95. }
  96.  
  97. gsl_dht *
  98. gsl_dht_new (size_t size, double nu, double xmax)
  99. {
  100.   int status;
  101.  
  102.   gsl_dht * dht = gsl_dht_alloc (size);
  103.  
  104.   if (dht == 0)
  105.     return 0;
  106.  
  107.   status = gsl_dht_init(dht, nu, xmax);
  108.   
  109.   if (status)
  110.     return 0;
  111.  
  112.   return dht;
  113. }
  114.  
  115. int
  116. gsl_dht_init(gsl_dht * t, double nu, double xmax)
  117. {
  118.   if(xmax <= 0.0) {
  119.     GSL_ERROR ("xmax is not positive", GSL_EDOM);
  120.   } else if(nu < 0.0) {
  121.     GSL_ERROR ("nu is negative", GSL_EDOM);
  122.   }
  123.   else {
  124.     size_t n, m;
  125.     int stat_bz = GSL_SUCCESS;
  126.     int stat_J  = 0;
  127.     double jN;
  128.  
  129.     if(nu != t->nu) {
  130.       /* Recalculate Bessel zeros if necessary. */
  131.       t->nu = nu;
  132.       stat_bz = dht_bessel_zeros(t);
  133.     }
  134.  
  135.     jN = t->j[t->size+1];
  136.  
  137.     t->xmax = xmax;
  138.     t->kmax = jN / xmax;
  139.  
  140.     t->J2[0] = 0.0;
  141.     for(m=1; m<t->size+1; m++) {
  142.       gsl_sf_result J;
  143.       stat_J += gsl_sf_bessel_Jnu_e(nu + 1.0, t->j[m], &J);
  144.       t->J2[m] = J.val * J.val;
  145.     }
  146.  
  147.     /* J_nu(j[n] j[m] / j[N]) = Jjj[n(n-1)/2 + m - 1], 1 <= n,m <= size
  148.      */
  149.     for(n=1; n<t->size+1; n++) {
  150.       for(m=1; m<=n; m++) {
  151.         double arg = t->j[n] * t->j[m] / jN;
  152.         gsl_sf_result J;
  153.         stat_J += gsl_sf_bessel_Jnu_e(nu, arg, &J);
  154.         t->Jjj[n*(n-1)/2 + m - 1] = J.val;
  155.       }
  156.     }
  157.  
  158.     if(stat_J != 0) {
  159.       GSL_ERROR("error computing bessel function", GSL_EFAILED);
  160.     }
  161.     else {
  162.       return stat_bz;
  163.     }
  164.   }
  165. }
  166.  
  167.  
  168. double gsl_dht_x_sample(const gsl_dht * t, int n)
  169. {
  170.   return t->j[n+1]/t->j[t->size+1] * t->xmax;
  171. }
  172.  
  173.  
  174. double gsl_dht_k_sample(const gsl_dht * t, int n)
  175. {
  176.   return t->j[n+1] / t->xmax;
  177. }
  178.  
  179.  
  180. void gsl_dht_free(gsl_dht * t)
  181. {
  182.   free(t->J2);
  183.   free(t->Jjj);
  184.   free(t->j);
  185.   free(t);
  186. }
  187.  
  188.  
  189. int
  190. gsl_dht_apply(const gsl_dht * t, double * f_in, double * f_out)
  191. {
  192.   const double jN = t->j[t->size + 1];
  193.   const double r  = t->xmax / jN;
  194.   size_t m;
  195.   size_t i;
  196.  
  197.   for(m=0; m<t->size; m++) {
  198.     double sum = 0.0;
  199.     double Y;
  200.     for(i=0; i<t->size; i++) {
  201.       /* Need to find max and min so that we
  202.        * address the symmetric Jjj matrix properly.
  203.        * FIXME: we can presumably optimize this
  204.        * by just running over the elements of Jjj
  205.        * in a deterministic manner.
  206.        */
  207.       size_t m_local; 
  208.       size_t n_local;
  209.       if(i < m) {
  210.         m_local = i;
  211.     n_local = m;
  212.       }
  213.       else {
  214.         m_local = m;
  215.     n_local = i;
  216.       }
  217.       Y = t->Jjj[n_local*(n_local+1)/2 + m_local] / t->J2[i+1];
  218.       sum += Y * f_in[i];
  219.     }
  220.     f_out[m] = sum * 2.0 * r*r;
  221.   }
  222.  
  223.   return GSL_SUCCESS;
  224. }
  225.  
  226.